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Abstract 

Factorization of polynomials is one of the foundations of symbolic com- 
putation. Its applications arise in numerous branches of mathematics and 
other sciences. However, the present advanced programming languages 
such as C++ and JH — h, do not support symbolic computation directly. 
Hence, it leads to difficulties in applying factorization in engineering fields. 
In this paper, we present an algorithm which use numerical method to ob- 
tain exact factors of a bivariate polynomial with rational coefficients. Our 
method can be directly implemented in efficient programming language 
such C++ together with the GNU Multiple-Precision Library. In addi- 
tion, the numerical computation part often only requires double precision 
and is easily parallelizable. 

Key words: Factorization of multivariate polynomials, Minimal poly- 
nomial, Interpolation methods, Numerical Continuation. 



1 Introduction 



Polynomial factorization plays a significant role in many problems including the 
simplification, primary decomposition, factorizcd Grobncr basis, solving poly- 
nomial equations and some engineering applications, etc. It has been studied 
for a long time and some high efficient algorithms have been proposed. There 
are two types of factorization approaches. One is the traditional polynomial 
factorization for exact input relying on symbolic computation, and the other is 
approximate polynomial factorization for inexact input. 
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t Corresponding author: dr.wenyuanwu@gmail.com 
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The traditional polynomial factorization methods follow Zassenhaus' approach 
[25] [26]. First, Multivariate polynomial factorization is reduced to bivariate fac- 
torization due to Bertini's theorem and Hensel lifting [5] [7] • Then one of the two 
remaining variables is specialized at random. The resulting univariate polyno- 
mial is factored and its factors are lifted up to a high enough precision. At last, 
the lifted factors are rccombined to get the factors of the original polynomial. 

Approximate factorization is a natural extension of conventional polynomial 
factorization. It adapts factorization problem to linear algebra first, then applies 
numerical methods to obtain an approximate factorization in complex which is 
the exact absolute factorization of a nearby problem. In 1985, Kaltofcn pre- 
sented an algorithm for computing the absolute irreducible factorization by 
floating point arithmetic |13j . Historically, the concept of approximate fac- 
torization appeared first in a paper on control theory [TS]. The algorithm is as 
follows: 1) represent the two factors G and H of the polynomials F with un- 
known coefficients by fixing their terms, 2) determine the numerical coefficients 
so as to minimize \\F — GH\\. Huang et al [8]. pursued this approach, but it 
seems to be rarely successful, unless G or H is a polynomial of several terms. 
In 1991, Sasaki et al. proposed a modern algorithm [19], which use power-series 
roots to find approximate factors. This algorithm is successful for polynomials 
of small degrees. Subsequently, Sasaki et al. presented another algorithm [20] 
which utilizes zero-sum relations. The zero-sum relations are quite effective for 
determining approximate factors. However, computation based on zero-sum re- 
lations is practically very time-consuming. In [21] . Sasaki presented an effective 
method to get as many zero-sum relations as possible by matrix operations so 
that approximate factorization algorithm is improved. Meanwhile, Corless et al 
proposed an algorithm for factoring bivariate approximate polynomial based on 
the idea of decomposition of affine variety in [4]. However, it is not so efficient 
to generalize the algorithm to multivariate case. A major breakthrough is due 
to Kaltofen et al. [12] [14] who extended Gao's work pj] from symbolics to 
numerics based on Ruppert matrix and Singular Value Decomposition. 

Symbolic factorization has been implemented in many Computer Algebra 
System. However, it is difficult to implemented directly in Programming Lan- 
guage such as C++ and J++, because most of Programming Language stan- 
dards do not support symbolic basic operators, and the compilers do not im- 
plement the symbolic computation, on which symbolic factorization is based. 
It restricts exact factorization from being applied in many engineering fields. 
Compared with symbolic factorization, approximate factorization can be imple- 
mented more easily in the popular programming languages. However it only 
gives approximate results even the input is exact. In this paper, we propose an 
almost completely numerical algorithm, which is not only implemented directly 
in the programming languages, but also achieve exact results. 

Except classic symbolic methods, some approaches have been proposed to 
obtain exact output by approximation [15] [ 27 ) . The idea of obtaining exact poly- 
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nomial factorization is from the connect between an approximate root of a given 
polynomial and its minimal polynomial in Q. Certainly, the minimal polyno- 
mial is a factor of the given input. Based on lattice basis reduced algorithm LLL 
and Integer Relation algorithm PSLQ of a vectors respectively, there are two 
algorithms for finding exact minimal polynomial of an algebraic number from its 
approximation. One is a numerical algorithm |15[ [5] for factorization of a uni- 
variate polynomial was provided by Transcendental Evaluation and high-degree 
evaluation, and the other for factorization of bivariate polynomial is based on 
LLL0I3]. But they are not efficient. 

In this paper, relying on LLL algorithm, we present an almost-completely 
numerical method for exact factoring polynomial with rational coefficient in 
Q. First, we choose a sample point in Q™ -1 at random. After specialization 
(i.e. substitution) at the point, the roots of the resulting univariate polynomial 
can be found very efficiently up to arbitrarily high accuracy. Then applying 
minimal polynomial algorithm to these roots yields an exact factorization of 
the univariate polynomial in Q. Next we shall move the sample point in "good 
direction" to generate enough number of points by using numerical continuation. 
Especially, for the rest sample points, the corresponding exact factorization can 
be found by using the same combination of roots as found in the first step. 
And these roots give more univariate polynomials for next step. Finally, the 
multivariate factorization can be obtained by interpolation. 

The paper is organized as follows. Section 2 gives a brief introduction of 
the preparation knowledge. Minimal polynomial algorithm will be discussed in 
Section 3. Then we present our method in Section 4,5 and 6. 



2 Preparation 

In this section, we briefly introduce the background knowledge and related topics 
to our article. 



2.1 Homotopy Continuation Methods 

Homotopy continuation methods play a fundamental role in Numerical Alge- 
braic Geometry and provide an efficient and stable way to compute all isolated 
roots of polynomial systems. These methods have been implemented in many 
software packages e.g. Hom4PS [IS], Bertini [T], PHCpack [2"4"] . 

The basic idea is to embed the target system into a family of systems con- 
tinuously depending on parameters. Then each point in the parameter space 
corresponds to a set of solutions. Suppose we know the solutions at a point. 
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Then we can track the solutions from this starting point to the point represent- 
ing the target system we want to solve. 



First let us look at the simplest case: a univariate polynomial f(z) with 
degree d. We know that f(z) has d roots in C (counting multiplicities). Of 
course we can embed f(z) into the family adz d + ad-iz^ 1 + • • • + do, where 
the di are parameters. Now choose a start point corresponding to z d — 1 in this 
parameter space, whose roots are 



_ 2kn y 

4>~U — C- 



* fc -c fc = 0,l,...,d-l (1) 

Then we use a real straight line in the parameter space to connect z d — 1 with 
/(*): 

H(z,t):=tf(z) + (l-t)(z d -l). (2) 
This form is a subclass of the family depending on only one real parameter 

ie [0,1]. 

When t = we have the start system H(z, 0) — z d — 1 and when t = 1 we 
have our target system H(z, 1) = /(z). An important question is to show how to 
track individual solutions as t changes from to 1. Let us look at the tracking of 
the solution z% (the fc-th root of f(z)). When t changes from to 1, it describes 
a curve, which is function of t, denoted by zy. = z^(t). So H(zk(t), t) = for all 
t G [0, 1]. Consequently, we have 

dH(z k (t),t) = dH( Zl t)dz k (t) dH(z,t) 

dt dz dt dt ' U 

This problem is reduced to an ode for the unknown function Zfc(t) together 
with an algebraic constraint H(zk(t),t) = 0. The initial condition is the start 
solution Zfc(0) = z\ and Zfc(l) is a solution of our target problem f(z) = 0. 



Remark 1 In the book J^, Blum, Smale et al. show that on average an approxi- 
mate root of a generic polynomial system can be found in polynomial time. Also 
application of the polynomial cost method for numerically solving differential 
algebraic equations \10jj gives polynomial cost method for solving homotopies. 



But there is a prerequisite for the continuous tracking: 9H q^ ^ along the 
curve z = z^{t). If the equations z — Zi-(t) = and tf'(z) + d(l — t)z d ~ x = 
have intersection at some point (t, zj,(t)), then we cannot continue the tracking. 
There is way to avoid this singular case, called the "gamma trick" that was first 
introduced in |17j . We know two complex curves almost always have intersec- 
tions at complex points, but here t must be real. So if we introduce a random 
complex transformation to the second curve, the intersection points will become 
complex points and such a singularity will not appear when t 6 [0, 1). Let us 
introduce a random angle 9 £ [— w, tt] and modify the homotopy ^) to 

H(z,t):=tf(z) + e ie (l-t)(z d -l). (4) 
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It is easy to show that the k-th starting solution is still z° in ((T|) and that Zfc(l) 
is still a root of f{z). 

2.2 Genericity and Probability One 

In an idealized model where paths arc tracked exactly and the random angle 
can be generated to infinite precision, the homotopy (jl]) can be proved to suc- 
ceed "with probability one". To clarify this statement, it is necessary to use a 
fundamental concept in algebraic geometry: genericity. 

Definition 1 (Generic) Let X be an irreducible algebraic variety. We say a 
property P holds generically on X , if the set of points of X that do not satisfy 
P are contained in a proper subvariety Y of X . The points in X\Y are called 
generic points. 

The set X\Y is called a Zariski open set of X. Roughly speaking, if Y is 
a proper subvariety of an irreducible variety X and p is a random point on X 
with uniform probability distribution, then the probability that p ^ Y is one. 
So we can consider a random point as a generic point on X without a precise 
description of Y. Many of the desirable behaviors of homotopy continuation 
methods rely on this fact. 

2.3 Coefficient-Parameter Homotopy 

There are several versions of the Coefficient-Parameter theorem in [22]. Here 
we only state the basic one. 

Theorem 1 Let F(z;q) = {fi(z; q), f n (z; q)} be a polynomial system in n 
variables z and m parameters q. Let Af(q) denote the number of nonsingular 
solutions as a function of q: 



1. There exist N, such that N(q) < N for any q € C m . Also {q G C m : 
J\f{q) = N} is a Zariski open set of C m . The exceptional set Y = {q : 
J\T(q) < N} is an affine variety contained in a variety with dimension 
m — 1. 




(5) 



Then, 
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2. The homotopy F(z; <j>(t)) = with 4>{t) : [0, 1) ->■ C m \y has N continuous 
non-singular solution paths z(t). 

3. When t — > 1~ , the limit of Zk(t), k = 1, ...,N includes all the non-singular 
roots of F(z; (f)(1)). 

An important question is how to choose a homotopy path <p(t) which can 
avoid the exceptional set Y. The following lemma [22] gives an easy way to 
address this problem. 

Lemma 1 Fix a point q and a proper algebraic set Y in C m . For a generic point 
p £ C m , the one-real- dimensional open line segment <p(t) := (1 — t) p + 1 q, t E 
[0, 1) is contained in C m \Y . 

2.4 Reductions 

Before factorization of a given polynomial, we shall first apply certain reduc- 
tions to the input to obtain a square-free polynomial over Q, which can remove 
multiplicities and ease the computation of the roots. Also we can assume each 
factor involves all the variables and has more than one term. Otherwise, we can 
compute the GCD to reduce the problem. For example, let F = f(x,y)g(y). 
Then F x = f x g and gcd(F, F x ) — g which gives us the factor g(y). 

By the Hilbert Irreducibility Theorem, we can further reduce the problem to 
univariate case by generic (random) specialization of one variable to a rational 
number. More precisely, if f(x,y) is irreducible in then for a generic 

rational number yo, f{x,yo) is also irreducible in the ring Q[x]. It means that 
the factorization is commutable with generic specialization. 

For univariate polynomial, there are symbolic methods to preform exact fac- 
torization in Q. Here we are more interested in numerical methods, namely 
from approximate roots to exact factors. 



3 Minimal Polynomial by Approximation 

There are two methods to compute the minimal polynomial of an algebraic 
number from its approximation. One is based on the LLL algorithm of the basis 
reduction [15], and another is based on PSLQ[5]. The later one is more efficient 
than the former one. However, it can only compute the minimal polynomial of 
a real algebraic number while the former one can find minimal polynomial of 
a complex algebraic number. Hence, we introduce the former algorithm which 
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is more suitable for this paper here. We refer the reader to the paper [T5] for 
more details. 

Let p(x) = Yll=o Pi x * be a polynomial. The length \p\ of p(x) is defined 
as the Euclidean norm of the vector (po,pi,--- ,p n ), and the height \p\oo as 
the Loo-norm of the vector (pg,pi, ■ • ■ ,Pn)- The degree and height of an alge- 
braic number are defined as the degree and height, respectively, of its minimal 
polynomial. 

Suppose that we have upper bound d and H on the degree and height re- 
spectively of an algebraic number with |a| < 1, and a complex rational number 
a approximating a such that \a\ < 1 and \a — d\ < 2~ s /(4d), where s is the 
smallest positive integer such that 

2 s > 2 d2 / 2 {d+l) {M+A)/2 H 2d 



Algorithm 1 [miniPoly] 

For n= 1, 2, • • • ,d in succession, do the following steps 



Step 1: construct 



( 1 

1 

1 

V 



2 s • Re(d ) 2 s ■ Im{a ) \ 

2 s -i?e(ai) 2 S -Im(ai) 

2 s • Re(a 2 ) 2 s ■ Im(a 2 ) 

1 2 s • Re(a n ) 2 s ■ Im(a n ) j 



(0) 



where Re(a) and Im(a) stand for the real part and imaginary part, 
respectively, of complex a, «o = 1 an d \di — a l \ < 2~ s ~ 1 / 2 for i = 
1, 2, ■ • • , n. Note 6>i can be computed by rounding the powers of a to 
s bits after the binary points. 

Step 2: Denote by bi the row i + 1 of the matrix in (0J). Apply the 
basic reduction algorithm to lattice L s = (bo, bi, • ■ ■ ,b n ), and obtain 
the reduced basis of the lattice. 

Step 3: If the first basis vector v = (vq, v\, ■ ■ ■ ,v n , v n +i, v n +2) in the 
reduced basis satisfies \v\ < 2 d ' 2 (d + 1)H , then return polynomial 
v i x ) = S™=o Vixl as rninimal polynomial of algebraic number a. 



Note: It is no major restriction to consider a with |a| < 1 only. In fact, 
if |a| > 1 satisfies the polynomial h(x) = 52j=o hiX*, then 1/a satisfies the 
polynomial Yli=o ^d-iX 1 ■ Therefore, if X^=o ^ l d-%x % is computed, the h(x) is 
obtained. Furthermore, an £-approximation a to a with |a| > 1 easily yields a 
3£-approximation /? to /3 = 1/a. This can be easily verified. 
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The following theorem shows the computation amount of calculating the 
minimal polynomial of an algebraic number [15j: 

Theorem 2 Let a be an algebraic number, and let d and H be upper bounds on 
the degree and height, respectively, of a. Suppose that we are given an approx- 
imation a to a such that \a — a\ < 2~ s /(12d), where s is the smallest positive 
integer such that 

2 s > 2 d2/2 (d+l) i3d+4)/2 H 2d 

Then the minimal polynomial of a can be determined in O(no ■ d 4 (d + logiJ)) 
arithmetic operations on integers having 0(cP • (<i + logiJ)) binary bits, where 
no is the degree of a. 



4 Finding More Polynomials by Continuation 

In the previous stage, we have the factors after specialization , which are univari- 
ate polynomials. To construct the factor of two variables by using interpolation, 
we need more information, i.e. specializations at more points. The main tool is 
the homotopy continuation method. 



4.1 Applying numerical continuation to factorization 

Suppose an input polynomial F(x, y) is reducible. Geometrically, if C denotes 
the zero set of / i.e. the union of many curves in C 2 , removing the singular locus 
of C from each curve Cj, the regular sets Si are connected in C 2 . Moreover, the 
singular locus has lower dimension, consequently it is a set of isolated points. 

Suppose f{x,y) is an irreducible factor of F in Q. Let yo,yi be random 
rational numbers. By the Hilbert Irreducibility Theorem the univariate polyno- 
mials fo = f(x,yo) and f\ = f(x,y\) are irreducible as well. Suppose we know 
the roots of fo. Then we can choose a path to connect yo and y\ avoiding the 
singular locus which has measure zero. By the Coefficient-Parameter Theorem, 
all the roots of f\ can be obtained by the following homotopy continuations: 

fix, y) = 

(l-t)(y-yo)+t(y-y 1 ) 1 = [) 

Moreover any generic complex number 7 implies that the homotopy path can 
avoid the singular locus by Lemma [T] when we track the path. 
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4.2 Control of the precision 



Let {xi, ..,x m } be the exact roots of fx and g be the primitive polynomial of 
fx. Then 

m 

g = a\\[x-Xi) £Z[i], (8) 

i=l 

for some integer number a. 

Note that we only have the approximate roots {xi, ..,x m }. 



Proposition 1 Letp = YiiLi( x ~ x i) o,ndp = YiiLx ( x ~ x i)- Let 5 = max, = i i .. im {|xj 
Xi\} and r = max, = i i .. im {|xi|}. If S is sufficiently small. Then 

||p-p||oo<( max {r'-'r^^m + l) S (9) 
V i=l,..,m \ l — 1 / / 



Proof. Let Xj = Xj + 5j. Thus, |<5j| < S. The left hand side ||p — p||oo = 
\\UiLi( x - x i+ 5 i)-IliLi(^-Xi)\\ = WHl'Li^^j^ ~ x j) S j\\+ °( S )- An upper 
bound of the coefficients of rii^ :) ( a:: — x j) with respect to x m ~ % is {f^Zi) 7 " 11 - 
Hence, ||p — p|[oo < max i= i > ... >m ("T! 1 )^ -1 • md + 6 □ 

Now let us consider how to find a. Suppose the input polynomial is F(x,y) 
and / is a factor of F. The primitive polynomial of /(x, y±), which is g, must 
be a factor of the primitive polynomial of F(x, yi). Thus, the leading coefficient 
of g must be a factor of the leading coefficient of the primitive polynomial of 
F(x, yi). Let a be the leading coefficient of the primitive polynomial of F(x, y\). 
Then let p = ctY[iLi( x — x i) ^ Z[x]. Note that itself may not be primitive, but 
its primitive polynomial is g. 

Let M = mnx, ,.. ,„{/ r 4 " 1 ^ 1 )} + 1. p = aUZi( x ~ x i)- Thus - if 6 < 
then ||p — p| |oo < 0.5. It means that we can round each coefficient of p to the 
nearest integer to obtain exact polynomial p which gives g. 



4.3 Detecting the degrees of factors 

After specialization at y = yo, we obtain the information about the number 
of factors and the degree of each factor with respect to x. The degrees with 
respect to y of factors provide the bound of the number of interpolation nodes. 
Certainly, we can use the degree of the input deg y (F) as the bound. However, 
the degrees with respect to y of factors are usually much less than deg y (F), 
especially when there are many factors. Therefore, for high efficiency, it is 
better to know the degree with respect to y of each factor. Now we will apply 
numerical continuation to detect such degree information. 
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We define the notation of 2-tuple degree to be 



deg(/) = [de g;c (/),deg y (/)]. 



Suppose deg(/) = [m, n] and / has r factors. Applying an approach of uni- 
variate polynomial solving to f(x,yo) and f{xo,y) yields points on the curve 
A = {(xi,y ), (x2,yo), • (x m ,yo)} and B = {(x , j/i), (x , y 2 ), (x , y„)} re- 
spectively. In addition, we also know the decomposition of two points sets in r 
groups with cardinalities {a%, ...,a r } and {61, b r }. Moreover J^a^ = m and 
2 h = n. 

Choose one point from each group of the first set A. Starting from these 
points, we track the homotopy path 



Because of the genericity of the choice of yo, xq and 7, the path avoids the 
singular locus. When t = 1, the endpoint must belong to the second set B. For 
example if the starting point of the first group of A and its end point belongs to 
the zth group of B. Then we know the first factor has degree [a\, hi]. Similarly, 
the degrees of other factors can be detected in the same way. 



5 Interpolation 

Polynomial interpolation is a classical numerical method. It is studied very 
well for univariate polynomials in numerical computation. Polynomial inter- 
polation problem is to determine a polynomial f(x) € F[x] with degree not 
greater than n £ N for a given pairs {(xi, fi), i = 0, • • ■ , n} satisfying f{x{) = /; 
for i = 0, • • • , n, where F is a field and Xj, fi G F. In general, there are four 
types of polynomial interpolation method: Lagrange Interpolation, Neville's 
Interpolation, Newton's Interpolation and Hermitc Interpolation. Lagrange in- 
terpolation and Newton's Interpolation formula are suited for obtaining inter- 
polation polynomial for a given set {(xj,/j),« = 0, • ■ • Neville's interpo- 
lation method aims at determining the value of the interpolating polynomial 
at some point. If the interpolating problem prescribes at each interpolation 
point {xi,i — 0, ■ • • , n} not only the value but also the derivatives of desired 
polynomial, then the Hermite formula is preferred. 

Different from the traditional interpolation problem above, our problem is to 
construct a bivariate polynomial from a sequence of univariate polynomials at 
chosen nodes. It is important to point out that the univariate polynomials are 
constructed by roots, which may not be equal to the polynomials by substitu- 
tions. But the only difference for each polynomial is just a scaling constant. 




(10) 
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More precisely, in this paper, we aim to solve a special polynomial interpola- 
tion problem: given a set of nodes and square free polynomials { (yi £ F, fi (x) € 
F[x}),i — 0, ■ • • , k}, compute a square free polynomial f(x, y) € F[x, y] of degree 
with respect to x not greater than n, where F is a field, such that f(x,yi) and 
fi{x) have the same roots. 

5.1 Illustrative examples 

EXAMPLE 5.1 Let f = x 2 + y 2 — 1. Since its degree with respect to y is 
two, we need three interpolation nodes which are y = —1/2,0,1/2. Suppose 
we know the roots at each node, then the interpolating polynomials are {/o = 
x 2 — 3/4, fi = x 2 — l,/2 = x 2 — 3/4}. To construct original polynomial f, we 
can use Lagrange method. 

Let t\ = ^2^)7-1/2- 1/2) ~ ~ y~ Similarly, £2 — — 4y 2 + 1 and £3 = 
2y 2 + y. It is easy to check that (x 2 - 3/4)^ + (x 2 - 1)£ 2 + (x 2 ~ 3/4)£ 3 = /. 

In the example above, the coefficient of / with respect to x 2 is a constant 1. 
Making the interpolating polynomials given by (jSJ) monic, we can construct / 
correctly by Lagrange basis. However, if the coefficient is nonconstant, i.e. a 
polynomial of y, then it is not straightforward to find /. The example below 
shows this problem. 

EXAMPLE 5.2 Let f = xy - 1. The nodes are y = 2,3. We know the roots 
are 1/2, 1/3 respectively at the nodes. Then the monic interpolating polynomials 
are {x — 1/2, x — 1/3}. If we still apply Lagrange basis £\ = —y + 3, £i~y — 2, 
it gives (x — l/2)(—y + 3) + (x — l/3)(y — 2) = x + 1/6 y — 5/6 which is totally 
different from the target polynomial xy — 1 . 

The basic reason is that the interpolating polynomials are not the polynomials 
after specialization s, and the only difference is certain scaling constants. To 
find these constants, we need more information. Now we use one more node: 
when y = 4, the monic interpolating polynomial is x — 1/4. By multiplying a 
scaling constant to / we can assume f(x, 4) = x — 1 /4, then there exist a, b such 
that f(x,2) = a(x— 1/2) and/(x,3) = b(x— 1/3). The corresponding Lagrange 
bases are i x = (y - 3)(y - 4)/2, £ 2 = -(y - 2)(y - 4), £ 3 = (y - 2){y - 3)/2. 
Then f = a (x- 1/2)4 + b(x - l/3)£ 2 + (x - 1/4)4. The coefficient of / with 
respect to y 2 must be zero. Consequently we have 

1/2 (x - 1/2) a + (x - 1/3)6 + 1/2 x- 1/8 = (11) 

which implies a linear system 

1/2 a -6 + 1/2 = 0,-1/4 a + 1/3 6- 1/8 = 
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The solution is a = 1/2, b = 3/4. Substituting them back to two nodes interpo- 
lation formula yields the polynomial we need, up to a constant 1 /4 



l/2(a; - l/2)(-y + 3) + 3/4(z - l/3)(y - 2) = (xy - l)/4 



5.2 Interpolation with indeterminates 



To extend the idea in example 15.21 we present a method to construct desired 
bivariate polynomial by using monic univariate interpolating polynomials. 

Suppose / is irreducible and its degrees with respect to x and y are m and n 
respectively. Consider x as the main variable, we can express this polynomial 
by / = J2iLo c i{y) xl i where c, are polynomials of y of degree less than or equal 
to n. We can consider each Cj as a vector in monomial basis. Suppose there 
are r linearly independent coefficients. If r = 1, then c,(y) = aiCo(y) for some 
constant a* and / = (Y]^L cnx 1 ) ■ co(y). It contradicts the assumption that / is 
irreducible. Hence, r > 2. 

Now we consider how to construct / by using the interpolating polynomi- 
als {fo(x),fi(x), ...,f k (x)} at k + 1 nodes {yo,yi, — ,J/k} respectively chosen at 
random. 

Let C be a (k + 1) x (m + 1) matrix [co,...,c m ] where is the column 
vector in monomial basis {y k ,y k ~ 1 , 1} of the polynomial Cj. Let V be the 



Vandermonde matrix 



Let A be a (fc + 1) x (m + 1) 



matrix where is the coefficient of the ith interpolating polynomial with 
respect to x 3 . To make the solution unique, we may fix f(x, yo) = /o and suppose 

Ai 

f(x, yi) = Xifi and A 4 7^ for i = 1, .., k. Let A 



V 



A fc / 



Therefore, 



V ■ C = A ■ A. 



(12) 



Since {yi} are distinct, the Vandermonde matrix has inverse and consequently 
C = V^ 1 ■ A • A. By our assumption, the degree with respect to y is n. It means 
that the first k — n rows of C must be zero. The zero at the ith row and jth 
column corresponds an equation. Thus, it leads to a linear system 



Row{V-\i) ■ A-Col{A,j) =0, 



(13) 
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for 1 < i < k — n and 1 < j < m + 1 with k unknowns. 

Only r linearly independent columns in A, so there are {k — n) r equations 
and k unknowns. The existence of the solution is due to the origination of the 
interpolating polynomials f(x, j/j) = Xifi for i = 1, k. The linear system has 
unique solution implies that (k — n)r > k. Thus, k > rn/(r — 1). Let u be the 
smallest integer greater than or equal to , namely 

" = r^i- (14) 



Thus, to determine the scaling constants {A^}, we need at least fj, more in- 
terpolation nodes. 

To find an upper bound for the number of nodes, let us consider / as 
a monic polynomial with rational function coefficients. All the coefficients 
{c m ,...,co} can be uniquely determined by rational function interpolation of 
X m + c m -i/c m x m ~ 1 + ■ ■ ■ + co/c m at 2n + 1 nodes. Therefore, it requires 2n 
nodes except the initial one. Thus, we have [i < k < 2n. 

But this upper bound is often overestimated, and for some special case the 
polynomial / can be constructed by using less nodes. 



Proposition 2 Let f be a polynomial in Q[x,y] and deg(/) = [m, n]. Suppose 
m > n and f has n + 1 linearly independent coefficients. Then f can be uniquely 
determined by n + 2 monic interpolating polynomials {fo{x), f n +i(x)} up to 
a scaling constant. 



Proof. Suppose the first n + 1 columns of A are linearly independent. By 
Equation (|12|) . we construct n + 1 equations: Row(V^ _1 , 1) ■ A-Co\(A,j) = 0, for 
j = 1, n + 1. Let B be the transpose of the submatrix consisting of the first 
n + 1 columns of A and v = (vi, ...,u„+2) t be the transpose of Row(V -1 , 1). 
Thus, 



= B- 



( vi 



•v = B- 



\ 



V 2 



•(Ai, A n+2 )* 



\ K+2 ) \ l'n+2 j 

Because the first n+1 coefficients of / are linearly independent, the evaluations 
of them at n + 2 random points must be linearly independent. So the rank 
of B is n + 1. Here v can be expressed by explicit form of the Vandcrmondc 
matrix [23| which is a vector of polynomials of {yo> Un+i}- For generic choice 
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of {yo, •••jj/ra+i}) each Vi ^ 0. Hence, the rank of B 



V2 



\ V n+2 J 

is still n + 1 and its nullity equals one. We can choose any solution {A^} to 
construct / by Lagrange basis: J^ILo ^ifi^i- ^ 

Remark 2 In our algorithm, we compute {Xi} starting from fi more nodes 
( together with the initial node j/o )> an d we add incrementally more nodes if 
necessary. But interestingly, the experimental results show that the lower bound 
H is often enough. This fact deserves further study. 



Algorithm 2 [Interpolation] 

Input : a set of polynomials {fo(x), ...,fk(x)} C Z[x, y] 

a set of rational numbers {yo, ■■■,yk} 

an integer n the degree of f with respect to y 
Output: f £ 1[x,y\, such that f(x,yi) = fi(x). 

1. Let A be the matrix consisting of the coefficient row vectors of the input 
univariate polynomials. 

2. Let r = Rank( A) . If k < /i, then it needs more interpolation nodes. 

3. Solve the homogenous linear system h!3\) to obtain the scaling constants 
{Ai, .., Xk} 



4- If the solution is not unique, then it needs more interpolation nodes. 

5. Else f = J2 k i=0 X i f i £ i €Q[x,y] 

6. Return the primitive polynomial of f 



6 Combination of Tools 



Now we combine the tools introduced in previous sections to obtain a new 
factorization algorithm. A preliminary version of the algorithm is implemented 
in Maple. For the efficiency, it requires a more sophisticate version in CH — h, 
even parallel program. 
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6.1 Main steps of the algorithm 

Algorithm 3 [Factorization] 
F = Fac(f) 

Input : /, a primitive polynomial / £ Z[x, y] such that gcd(/, f x ) = 1 
Output: F, a set of primitive polynomials {/i, f r } C 1\x, y\, such that / 

1 . Apply a numerical solver to approximate the roots of f{x,yo) = and 
f(x ,y) = at generic points x 0l y a E Q. 

2. Apply miniPoly to roots above and decompose the solutions. And gener- 
ate the minimal polynomials for them and we have grouping information 
of roots. In this step, it needs Newton iteration to refine the roots up to 
desired accuracy. 

3. Apply homotopy (flU)) to obtain the degrees of each factor. 

4. For group i (corresponding to the factor /j), i = 1, r, use homotopy (J7J 
to generate its approximate roots at random rational numbers {yi, yk}- 

5. For each set of roots at yj, refine the roots to the accuracy given by Propo- 
sition [1] then make the product and construct the polynomial fi(x,yj). 

6. Call interpolate with the interpolating polynomials {fi(x, yo), fi{x, yk)} 
to construct fi(x,y). 

6.2 A simple example 

Let us consider a polynomial / = (xy — 2) (x 2 + y 2 — 1). First, we choose a 
sequence of random rational numbers {97/101, 1, 104/101, 123/101, 129/101, ...}. 
Substituting y = 97/101 into / yields Mignotte bound of the coefficients of 
factors 9170981 and the digits required to produce the minimal polynomial is 110 
by Theorem [2] Then compute the approximate roots of /(x, 97/101) up to 110 
digits accuracy. The miniPoly subroutine gives two groups of points: [[1, 2], [3]] 
and the corresponding minimal polynomials [—792 + 10201 x 2 , —202 + 97 a:]. By 
Hilbert Irreducibility theorem, there are two factors. On the other hand, fix the 
value of x and obtain the univariate polynomials [—202 + 97 y , —792 + 10201 y 2 } 
and [[3], [1,2]]. 

Starting from the first point of group one, the Homotopy (JTUJ) path ends at 
a point which satisfies y, -792 + 10201 y 2 . It implies that -792 + 10201 x 2 and 
y, -792 + 10201 y 2 arc from the same factor of degree [2, 2]. By Equation ([13]). 
we need /i = [^rr] = 4 more interpolating polynomials which are produced 
by Homotopy 0. Thus, there are five polynomials [-792 + 10201 x 2 , x 2 , 615 + 
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10201 x 2 , 4928 + 10201 x 2 , 6440 + 10201 x 2 ]. The scaling constants [Ai = 1, A 2 = 
10201, A3 = 1,A 4 = 1,A 5 = 1] are obtained by system (fT3|) . Consequently, the 
Lagrange interpolation formula gives the correct factor — 1 + x 2 + y 2 . 

Since the degree of the other factor is [1, 1], it needs fi = 2 more polynomi- 
als and they arc [—202 + 97 x,x — 2, —101 + 52 x\. The corresponding scaling 
constants are [Ai = 1/2, A2 = 101/2, A3 = 1] and the resulting factor is x y — 2. 



7 Conclusion 

A new numerical method to factorize bivariate polynomials exactly is presented 
in this article. We implemented our algorithm in Maple to verify the correct- 
ness. More importantly, the main components of our algorithm, miniPoly and 
numerical homotopy continuation can be implemented directly in CH — h or JH — h 
with existing multi-precision packages, e.g. GNU MP library. Furthermore, 
these two numerical components are naturally parallclizablc. Therefore, it gives 
an alternative way to exact factorization which can take the advantages of stan- 
dard programming languages and parallel computation techniques widely used 
by industries. 

In this article, we mainly focus on bivariate case. It is quite straightforward 
to extend to multivariate case. However, the number of the interpolation nodes 
grows exponentially as the increasing of the number of monomials. A more 
practical way to deal with such difficulty is to exploit the sparsity if the factors 
are sparse. It desires the further study in our future work. 
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